/*
 * Copyright (c) Hunan University.
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

#include <iostream>
#include "nlcpp.h"

using namespace nlcpp;
using namespace std;

void test_linear_equation_solver() {
    Matrix<double> A(3, 3);
    Vector<double> b(3);

    A << 10, 3, 1, 
        2, -10, 3, 
        1, 3, 10;
    b << 14, -5, 14;

    cout << "========Test Linear Equation Solvers (Case 1) ...========" << endl;
    
    auto jacobi_iteration_result = LinearEquationSolver::jacobi_iteration(A, b, 0.01);
    cout << "Result by jacobi iteration " << endl << jacobi_iteration_result << endl;

    auto gauss_seidel_iteration_result = LinearEquationSolver::gauss_seidel_iteration(A, b, 0.01);
    cout << "Result by gauss-seidel iteration " << endl << gauss_seidel_iteration_result << endl;

    auto gauss_elimination_result = LinearEquationSolver::gauss_elimination(A, b);
    cout << "Result by gauss elimination " << endl << gauss_elimination_result << endl;

    auto lu_decomposition_result = LinearEquationSolver::solve_lu_decomposition(A, b);
    cout << "Result by lu decomposition " << endl << lu_decomposition_result << endl;
    
    cout << endl;

    Matrix<double> A2(4, 4);
    Vector<double> b2(4);

    A2 << 6, 2, 1, -1,
        2, 4, 1, 0,
        1, 1, 4, -1,
        -1, 0, -1, 3;
    b2 << 6, -1, 5, -5;

    cout << "========Test Linear Equation Solvers (Case 2) ...========" << endl;

    jacobi_iteration_result = LinearEquationSolver::jacobi_iteration(A2, b2, 0.01);
    cout << "Result by jacobi iteration " << endl << jacobi_iteration_result << endl;

    gauss_seidel_iteration_result = LinearEquationSolver::gauss_seidel_iteration(A2, b2, 0.01);
    cout << "Result by gauss-seidel iteration " << endl << gauss_seidel_iteration_result << endl;

    gauss_elimination_result = LinearEquationSolver::gauss_elimination(A2, b2);
    cout << "Result by gauss elimination " << endl << gauss_elimination_result << endl;

    lu_decomposition_result = LinearEquationSolver::solve_lu_decomposition(A2, b2);
    cout << "Result by lu decomposition " << endl << lu_decomposition_result << endl;

    cout << endl;

    Matrix<double> A3(4, 4);
    Vector<double> b3(4);

    A3 << 1, 0, 0, 0,
        1.0/3.0, 1, 0, 0,
        1.0/6.0, 1.0/5.0, 1, 0,
        -1.0/6.0, 1.0/10.0, -9.0/37.0, 1;
    b3 << 6, -1, 5, 5;

    cout << "========Test Linear Equation Solvers (Case 3) ...========" << endl;

    jacobi_iteration_result = LinearEquationSolver::jacobi_iteration(A3, b3, 0.01);
    cout << "Result by jacobi iteration " << endl << jacobi_iteration_result << endl;

    gauss_seidel_iteration_result = LinearEquationSolver::gauss_seidel_iteration(A3, b3, 0.01);
    cout << "Result by gauss-seidel iteration " << endl << gauss_seidel_iteration_result << endl;

    gauss_elimination_result = LinearEquationSolver::gauss_elimination(A3, b3);
    cout << "Result by gauss elimination " << endl << gauss_elimination_result << endl;

    lu_decomposition_result = LinearEquationSolver::solve_lu_decomposition(A3, b3);
    cout << "Result by lu decomposition " << endl << lu_decomposition_result << endl;

    cout << endl;
}
